{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transport measurement data analysis\n",
    "This is an example notebook for the analysis class `IV_curve` of `qkit.analysis.IV_curve.py`. This handels transport measurment data (focussed of measurements of Josephson junctions in the current bias) taken with `qkit.measure.transport.transport.py` and provides methods to\n",
    "* load data files,\n",
    "* merge data files,\n",
    "* calculate numerical derivatives, such as differential resistances,\n",
    "* analyse voltage and current offsets,\n",
    "* correct ohmic resistances offsets arising due to lead resistivity in 2-wire measurements,\n",
    "* analyse the normal state resistance,\n",
    "* analyse critical currents and voltage jumps,\n",
    "* analyse switching current distributions.\n",
    "\n",
    "For error propagation the [`uncertainties`](https://github.com/lebigot/uncertainties) package is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from uncertainties import ufloat, umath, unumpy as unp\n",
    "from scipy import signal as sig\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qkit\n",
    "qkit.start()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qkit.analysis.IV_curve import IV_curve as IVC\n",
    "ivc = IVC()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load qkit transport measurement file\n",
    "Transport measurement data with a given `uuid` can be loaded using `ivc.load(uuid)`. Several elements are available, especially\n",
    "* data file `ivc.df`,\n",
    "* settings `ivc.settings`,\n",
    "* measurement object `ivc.mo`,\n",
    "* current values `ivc.I`,\n",
    "* voltage values `ivc.V`,\n",
    "* differential resistance values `ivc.dVdI`,\n",
    "* sweep list `ivc.sweeps`,\n",
    "* scan dimension (1D, 2D or 3D) `ivc.scan_dim`,\n",
    "* in case of 2D and 3D scans, x-parameter dataset `ivc.x_ds`, values `ivc.x_vec`, name `ivc.x_coordname`, unit `ivc.x_unit`, \n",
    "* in case of 3D scans, y-parameter dataset `ivc.y_ds`, values `ivc.y_vec`, name `ivc.y_coordname`, unit `ivc.y_unit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.load(uuid='XXXXXX')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge qkit transport measurement files\n",
    "Qkit transport measurement files can be merged depending on the scan dimension to one single new file by `ivc.merge()`.\n",
    "* 1D: all sweep data are stacked and views are merged.\n",
    "* 2D: values of x-parameter and its corresponding sweep data are merged in the order `order`.\n",
    "* 3D: values of x- and y-parameters and its corresponding sweep data are merged in the order `order`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.merge(uuids=('XXXXXX', 'YYYYYY'), order=(-1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential resistance\n",
    "The differential resistance $ \\frac{\\text{d}V}{\\text{d}I} $ is caclulated as numerical derivative using `ivc.get_dVdI()`. By default the Savitzky Golay filter is applied, but different methods can be used, e.g. a simple numerical gradient `ivc.get_dVdI(mode=np.gradient)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.get_dVdI()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Current and voltage offsets\n",
    "The current and voltage offsets can be calculated using `ivc.get_offsets()` or `ivc.get_offset()`.\n",
    "The branch where the y-values are nearly constant are evaluated. The average of all corresponding x-values is considered to be the x-offset and the average of the extreme y-values are considered as y-offset. These are by default the critical y-values $ y_c$, but can also be set to retrapping y-values $ y_r $ if `yr=True`\n",
    " * `ivc.get_offsets()` calculates x- and y-offset of every trace,\n",
    " * `ivc.get_offset()` calculates x- and y-offset of the whole set (this differs only for 2D or 3D scans).\n",
    " \n",
    "Note that reasonable initial values `offset` and `tol_offset` are sufficient to find the range where the y-values are nearly constant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.get_offsets(offset=0, tol_offset=20e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.get_offset(offset=0, tol_offset=20e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ohmic resistance offset\n",
    "The voltage values can be corrected by an ohmic resistance offset such as occur in 2wire measurements using `ivc.get_2wire_slope_correction()`. The two maxima in the differential resistivity $ \\frac{\\text{d}V}{\\text{d}I} $ are identified as critical and retrapping currents. This is done using `scipy.signal.find_peaks()` by default, but can set to custom peak finding algorithms `peak_finder`. The slope of the superconducting regime in between (which should ideally be infinity) is fitted using `numpy.linalg.qr()´ algorithm and subtracted from the raw data.\n",
    "\n",
    "Note that the arguments of the peak finding algorithms need to be set properly, e.g. `prominence` for `scipy.signal.find_peaks()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.get_2wire_slope_correction(prominence=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Critical and retrapping currents\n",
    "Critical currents can be determined on three different ways. The voltage jump at the critical and retrapping current can be found by\n",
    "* a voltage threshold value that is exceeded,\n",
    "* peaks in the differential resistance $ \\frac{\\text{d}V}{\\text{d}I} $,\n",
    "* peaks in the Gaussian smoothed derivative $ \\text{i}f\\exp\\left(-sf^2\\right) $ in the frequency domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Voltage threshold methode for critical and retrapping currents\n",
    "The critical and retrapping currents $ I_\\text{c} $, $ I_\\text{r} $ can be calculated by finding the currents that correspond to the voltages which exceed a certain threshold using `ivc.get_Ic_threshold()`. The branch where the voltage values are nearly constant are evaluated. Their maximal values of the up- and down-sweep are considered as critical currents $ I_c $ and retrapping current $ I_r $ (if `Ir=True`), respectively.\n",
    "\n",
    "Note that it works best, if `offset` is already determined via `get_offsets()` and that a reasonable initial value `tol_offset` is sufficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.get_Ic_threshold(Ir=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Peak detection in differential resistance methode for critical and retrapping currents\n",
    "The critical and retrapping currents $ I_\\text{c} $, $ I_\\text{r} $ can be calculated by detecting peaks in the differential resistance $ \\frac{\\text{d}V}{\\text{d}I} $ using `ivc.get_Ic_deriv()`. This is done using `scipy.signal.find_peaks()` by default, but can set to custom peak finding algorithms `peak_finder`.\n",
    "\n",
    "Note that the arguments of the peak finding algorithms need to be set properly, e.g. `prominence` for `scipy.signal.find_peaks()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_cs, I_rs, props = ivc.get_Ic_deriv(prominence=1, Ir=True)\n",
    "I_cs, I_rs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('all currents, where voltage jumps')\n",
    "if ivc._scan_dim == 1:\n",
    "    print(np.array(map(lambda p1D: p1D['I'], props)))\n",
    "elif ivc._scan_dim == 2:\n",
    "    print(np.array(map(lambda p2D: map(lambda p1D: p1D['I'], p2D), props)))\n",
    "elif ivc._scan_dim == 3:\n",
    "    print(np.array(map(lambda p3D: map(lambda p2D: map(lambda p1D: p1D['I'], p2D), p3D), props)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Peak detection in the Gaussian smoothed derivative methode for critical and retrapping currents\n",
    "The critical and retrapping currents $ I_\\text{c} $, $ I_\\text{r} $ can be calculated by detecting peaks in the Gaussian smoothed derivative $ \\left(\\text{i}f\\cdot\\text{e}^{-sf^2}\\right) $ in the frequency domain using `ivc.get_Ic_dft()`. This is done using `scipy.signal.find_peaks()` by default, but can set to custom peak finding algorithms `peak_finder`.\n",
    "\n",
    "Note that the smoothing factor `s` and the arguments of the peak finding algorithms need to be set properly, e.g. `prominence` for `scipy.signal.find_peaks()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_cs, I_rs, props = ivc.get_Ic_dft(Ir=True, prominence=1e-6)\n",
    "I_cs, I_rs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normal state resistance \n",
    "The normal state resistance $ R_\\text{n} $ of ohmic (overcritical) branch can be calculated using `ivc.get_Rn()`.\n",
    "The normal state resistance corresponds to the inverse linear slope of the normal conducting branch $ I=R_\\text{n}^{-1}U $ (`mode=0`) or the average value of $ \\frac{\\mathrm{d}V}{\\mathrm{d}I} $ of the normal conducting branch (`mode=1`). The ohmic range, in turn, is considered to range from the outermost tail of the peaks in the curvature $ \\frac{\\mathrm{d}^2V}{\\mathrm{d}I^2} $ to the start/end of the sweep and the resistance is calculated as mean of the differential resistance values `dVdI` within this range. This is done using `scipy.signal.savgol_filter(deriv=2)` and `scipy.signal.find_peaks()` by default, but can set to any second order derivative function `deriv_func` and peak finding algorithms `peak_finder`. For `scipy.signal.find_peaks()` the `prominence` parameter is set to $ 1\\,\\% $ of the absolute curvature value by default.\n",
    "\n",
    "Note that the arguments of the peak finding algorithms need to be set properly, e.g. `prominence` for `scipy.signal.find_peaks()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ivc.get_Rn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Switching current measurements\n",
    "Switching current distributions can be analyzed and plotted using the `ivc.scm.fit()` and `ivc.scm.plot()` of the `switching_current` subclass.\n",
    "\n",
    "The switching currents need to be determined beforehand, such as by `ivc.get_Ic_deriv()`. Their switching current distribution $ P_k \\mathrm{d}I = \\frac{n_k}{N\\Delta I}\\mathrm{d}I $ is normalized $ \\int\\limits_0^\\infty P(I^\\prime)\\mathrm{d}I^\\prime = 1 $ and get by `numpy.histogram()`.\n",
    "The escape rate reads $ \\Gamma(I_k) = \\frac{\\left|\\frac{\\mathrm{d}I}{\\mathrm{d}t}\\right|}{\\Delta I}\\ln\\left(\\frac{\\sum\\limits_{j\\geq k} P_j}{\\sum_\\limits{j\\geq k+1} P_j}\\right) $ and the normalized escape rate $ \\left[\\ln\\left(\\frac{\\omega_0}{2\\pi\\Gamma}\\right)\\right]^{2/3} $ is fitted versus $ \\gamma $ to $ f(\\bar{\\gamma}) = a\\cdot\\bar{\\gamma}+b $ where the root $ \\left[\\ln\\left(\\frac{\\omega_0}{2\\pi\\Gamma}\\right)\\right]^{2/3} = 0 $ yields the critical current $ I_\\text{c} = -\\frac{b}{a} $. Here, the sweep rate results as $ \\frac{\\mathrm{d}I}{\\mathrm{d}t} = \\delta I\\cdot\\frac{\\text{nplc}}{\\text{plc}} $, the centers of bins as moving average of the returned bin-edges using `np.convolve(edges, np.ones((2,))/2, mode='valid')` and the bin width as $ \\Delta I = \\frac{\\max(I_\\text{b})-\\min(I_\\text{b})}{N_\\text{bins}} $. \n",
    "\n",
    "### theoretical background\n",
    "* the probability distribution of switching currents is related to the escape rate $ \\Gamma(I) $ and the current ramp rate $ \\frac{\\mathrm{d}I}{\\mathrm{d}t} $ as\n",
    "\\begin{equation}\n",
    "P(I)\\mathrm{d}I = \\Gamma(I)\\left|\\frac{\\mathrm{d}I}{\\mathrm{d}t}\\right|^{-1} \\left(1-\\int\\limits_0^I P(I^\\prime)\\mathrm{d}I^\\prime\\right)\\mathrm{d}I\n",
    "\\end{equation}\n",
    "This integral equation can be solved explixitly for the switching-current distribution\n",
    "\\begin{align}\n",
    "P(I) &= \\Gamma(I)\\left|\\frac{\\mathrm{d}I}{\\mathrm{d}t}\\right|^{-1}\\exp\\left(-\\left|\\frac{\\mathrm{d}I}{\\mathrm{d}t}\\right|^{-1}\\int\\limits_0^I\\Gamma(I^\\prime)\\mathrm{d}I^\\prime\\right)\\\\\n",
    "P(I_k) &= \\Gamma(I_k)\\left|\\frac{\\mathrm{d}I_k}{\\mathrm{d}t}\\right|^{-1}\\exp\\left(-\\left|\\frac{\\mathrm{d}I_k}{\\mathrm{d}t}\\right|^{-1}\\Delta I\\sum\\limits_{j=0}^k\\Gamma(I_j)\\right)\n",
    "\\end{align}\n",
    "* Solving for the escape rate results in\n",
    "\\begin{align}\n",
    "\\Gamma(I) &= \\frac{\\left|\\frac{\\mathrm{d}I}{\\mathrm{d}t}\\right|}{\\Delta I}\\ln\\left(\\frac{\\int\\limits_I^\\infty P(I^\\prime)\\mathrm{d}I^\\prime}{\\int\\limits_{I+\\Delta I}^\\infty P(I^\\prime)\\mathrm{d}I^\\prime}\\right) \\\\\n",
    "\\Gamma(I_k) &= \\frac{\\left|\\frac{\\mathrm{d}I}{\\mathrm{d}t}\\right|}{\\Delta I}\\ln\\left(\\frac{\\sum\\limits_{j\\geq k} P_j}{\\sum_\\limits{j\\geq k+1} P_j}\\right)\n",
    "\\end{align}\n",
    "* The escape rate, in turn, is related to the attempt frequency $ \\frac{\\omega_0}{2\\pi} $ and the barrier height $ U_0 = 2E_\\mathrm{J}\\left(\\sqrt{1-\\gamma^2}-\\gamma\\arccos(\\gamma)\\right) \\approx E_\\mathrm{J}\\frac{4\\sqrt{2}}{3}\\left(1-\\gamma\\right)^\\frac{3}{2} $ and results in\n",
    "\\begin{align}\n",
    "\\Gamma_\\text{th} &= \\frac{\\omega_0}{2\\pi}\\exp\\left(\\frac{U_0}{k_\\text{B}T}\\right) \\\\\n",
    "&= \\frac{\\omega_0}{2\\pi}\\exp\\left(-\\frac{E_\\text{J}\\frac{4\\sqrt{2}}{3}(1-\\bar{\\gamma})^{3/2}}{k_\\text{B}T}\\right)\\\\\n",
    "\\left[\\ln\\left(\\frac{\\omega_0}{2\\pi\\Gamma}\\right)\\right]^{2/3} &= \\left(\\frac{E_\\text{J}}{k_\\text{B}T}\\frac{4\\sqrt{2}}{3}\\right)^{2/3}\\cdot(1-\\bar{\\gamma})\n",
    "\\end{align}\n",
    "* References\n",
    "  * [Fulton, T. A., and L. N. Dunkleberger. \"Lifetime of the zero-voltage state in Josephson tunnel junctions.\" Physical Review B 9.11 (1974): 4760.](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.9.4760)\n",
    "  * [Wallraff, Andreas. Fluxon dynamics in annular Josephson junctions: From Relativistic Strings to quantum particles. Lehrstuhl für Mikrocharakterisierung, Friedrich-Alexander-Universität, 2001.](https://www.osti.gov/etdeweb/biblio/20203466)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "props = ivc.get_Ic_deriv(prominence=1)\n",
    "I_0 = np.array(list(map(lambda p2D: list(map(lambda p1D: p1D['I'][0], p2D)), props)))\n",
    "ivc.scm.fit(I_0=I_0*1e6,\n",
    "            omega_0=1e9,\n",
    "            bins=30)\n",
    "ivc.scm.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
